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We present a nonperturbative computation of the equation of state of polarized, attractively in¬ 
teracting, nonrelativistic fermions in one spatial dimension at finite temperature. We show results 
for the density, spin magnetization, magnetic susceptibility, and Tan’s contact. We compare with 
the second-order virial expansion, a next-to-leading-order lattice perturbation theory calculation, 
and interpret our results in terms of pairing correlations. Our lattice Monte Carlo calculations im¬ 
plement an imaginary chemical potential difference to avoid the sign problem. The thermodynamic 
results on the imaginary side are analytically continued to obtain results on the real axis. We focus 
on an intermediate- to strong-coupling regime, and cover a wide range of temperatures and spin 
imbalances. 


I. INTRODUCTION 

In the last two decades, experimental studies with ul¬ 
tracold atomic gases have made consistent strides to¬ 
wards the increasingly clean and controlled characteriza¬ 
tion of strongly coupled matter, particularly in nonper¬ 
turbative regimes that are out of reach for conventional 
theory methods Q]. This represents a challenge to the 
theory side, which has been met in some cases but re¬ 
mains open in general. 

As is well known, spatial dimensionality plays a cru¬ 
cial role in these systems regardless of the strength of 
the interaction. Although it is also generally understood 
that interactions tend to dominate in lower dimensions 
and, conversely, mean-field descriptions become more re¬ 
liable in higher dimensions, they in general do not al¬ 
low for quantitative predictions. This can be viewed 
as a signal that fluctuation effects still play a promi¬ 
nent role, calling for more sophisticated theoretical tools. 
Both from the theory and experiment sides, considerable 
progress has been made in the study of three-dimensional 
(3D) systems in a variety of situations (e.g., in harmonic 
traps, homogeneous space, polarized, unpolarized, in the 
ground state, at finite temperature, etc.; see Ref. [2] for 
reviews), in particular with emphasis on the so-called 
BEC-BCS crossover and scale-invariant regimes [3j, as 
well as the Efimov effect [3]. In recent years, there has 
also been increasing activity in similar studies in 2D (see, 
e.g., Ref. [5] for a recent review), where the possibility of 
accessing directly the superfluid Berezinskii-Kosterlitz- 
Thouless transition ,B| has been a major drive. 

In this context, the motivation to study one¬ 
dimensional systems, in particular fermions, is manifold. 
Large classes of ID problems can be solved exactly at zero 
temperature via powerful techniques such as the Bethe 
ansatz, which has propelled a fair amount of work over 
the last few decades (see, e.g., Refs. El). However, a 
new wave of interest has been underway for a few years. 
This renewed activity stems in part from the realization 
of ID systems in the form of atomic gases in highly con¬ 
strained, quasi-ID optical traps, but it is also due to the 


advent of quantum-information concepts in condensed- 
matter physics 1 and their connection to quantum phase 
transitions (in particular topological ones) in low dimen¬ 
sions. 

Systems of spin-polarized, attractively interacting 
fermions are particularly appealing because of the po¬ 
tential occurrence of spontaneously broken translation in¬ 
variance in exotic superfluid phases nnj. The focus of the 
present work, however, is on the basic thermodynamic 
equations of state of spin-polarized fermions in ID, rather 
than on detecting potential exotic superfluid phases. We 
compute the density and the spin magnetization as well 
as Tan’s contact, which encodes the importance of high- 
momentum correlations in systems with short-range in¬ 
teractions mm ■ Our study is distinguished from previ¬ 
ous ones in that we make use of complex chemical poten¬ 
tials to overcome the so-called sign problem. The latter 
has been a major roadblock in lattice Monte Carlo calcu¬ 
lations of asymmetric systems (e.g., systems with mass 
or spin imbalance mm, or at finite quark density in 
the case of QCD m as explained below. As a proof 
of principle, we consider spin-1/2 fermions in ID with 
attractive contact interactions, but more general situa¬ 
tions including richer interactions and higher dimensions 
can be studied with the same methods. Note that we 
do not determine the exact nature of the ground state of 
the theory in our present study; as we explain in detail 
below, our approach based on complex chemical poten¬ 
tials is in fact unable to reach the ground state (and at 
fixed temperature is not able to reach arbitrary polariza¬ 
tions). Nevertheless, the ground state indirectly leaves its 
imprint in our computation of the basic thermodynamic 
equations of state. 

The Hamiltonian we analyze is that of the Gaudin- 
Yang model m, 

i i<j 

where the sums are over all particles. In our study, we 
will consider polarized systems, in the sense that they 
will have a nonvanishing average spin magnetization in 
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general. In the grand-canonical ensemble, the partition 
function of such a system is 


£(0^,0^) = Trex P ~ ~ M 4 .-N 4 J > ( 2 ) 


where /i s is the chemical potential for spin s =j\ i par- 
tides, N s is the corresponding particle number operator, 
and /3 is the inverse temperature. 

As mentioned above, asymmetric systems are chal¬ 
lenging for lattice Monte Carlo calculations due to the 
appearance of the sign problem. To circumvent this 
difficulty, we implement the approach put forward in 
Ref. [H]. Here we attempt this type of calculation for a 
nonrelativistic theory; a similar strategy has been applied 
by some of the present authors to ground-state calcula¬ 
tions in the mass-imbalanced case m. In the present 
case, we take the chemical potential for each fermion 
species to be complex, but such that one is the com¬ 
plex conjugate of the other: m = /zT. As a consequence, 
the fermion determinants in the Monte Carlo calculation 
(see, e.g., Ref. [HI [T7j) are also complex conjugates of 
one another, and the probability measure is thus non¬ 
negative. 

In this approach, the overall chemical potential fi = 
(//| + /i, )/2 is real, as usual, but the asymmetry param¬ 
eter h = (/r^. — g^)/2 is imaginary. For convenience, we 
define hi := Im h = —ih, with h\ being a real-valued 
number. The total density n = n* + n , in our study 
based on the grand-canonical ensemble is then still a real¬ 
valued number, while the so-called spin magnetization 
rri = ru — n^ is imaginary. These, as well as every output 
of the calculation, must be analytically continued to the 
real -h axis in order to obtain the physical results. While 
the analytic continuation procedure introduces some de¬ 
gree of arbitrariness in the final results (see below), it 
should be pointed out that the results on the imaginary 
side are fully nonperturbative and, in principle, exact, 
and certain aspects of the functional dependence with 
respect to h are known. Note, for instance, that the 
asymmetry h always enters in the calculations as a func¬ 
tion of f3h. Moreover, it can be shown that the results 
for imaginary asymmetries will be 27 t periodic in /?/q; see 
Ref. [13]. This is therefore a compact parameter and we 
will restrict it to the interval [— 7 r, 7 r]. The symmetry in 
the form of the partition function under spin exchange 
indicates that the physics it describes is independent of 
the sign of h 7 i.e., we expect our results to be either odd 
or even functions of h, depending on the observable. 

We would like to emphasize that although it is, in prin¬ 
ciple, exact, our present Monte Carlo approach to spin- 
imbalanced Fermi gases is not capable of studying the 
zero-temperature limit, but is limited to finite temper¬ 
atures \f3hi\ < 7 r being a direct consequence of the 27 t 
periodicity in j3hi. For the same reason, one may say 
that at fixed temperature (i.e. /3) not all polarizations 
are achievable, as the above constraint limits the range 
of hi. In this sense, our present approach is complemen¬ 
tary to the Bethe ansatz 0 i which allows for an ex¬ 


act solution of one-dimensional Fermi gases in the zero- 
temperature limit. Still, our approach enables us to com¬ 
pute the finite-temperature equation of state in a certain 
parameter range which can then potentially be compared 
to experiments, see, e.g., Ref. m- The present, however, 
should rather be regarded as a proof-of-principle work 
that serves as an intermediate step towards calculations 
in higher dimensions. 

It should also be pointed out that we do not expect 
phase transitions to be present in this ID system as a 
function of temperature due to the absence of sponta¬ 
neous symmetry breaking m■ This fact should allow for 
a more reliable analytic continuation from /?hi to the real 
side. Conversely, the appearance of accumulation points 
of zeros of the partition function (the so-called Yang-Lee 
zeros) in higher dimensions may complicate the analytic 
continuation in those cases. 


II. SCALES AND COMPUTATIONAL 
TECHNIQUE 

A. General setup 

The problem of fermions with a contact interaction is 
ultraviolet-finite in one spatial dimension. Therefore, the 
bare coupling has a physical meaning: in the continuum 
limit, g = 2/ao, where ao is the scattering length for the 
symmetric channel (see e.g. [ 20 ] j: accordingly, g will be 
reported in units of y/]3. Note that the thermal de Broglie 
wavelength is \ T = 7]3. 

To characterize the thermodynamics of polarized in¬ 
teracting fermions in one dimension we compute three 
quantities, namely the density n, the spin magnetization 
m and the contact C, as functions of the inverse temper¬ 
ature /3, the average chemical potential /x, the chemical 
potential difference h, and the coupling g. To make all 
of these quantities dimensionless, we utilize the noninter¬ 
acting unpolarized density n 0 as the scale for n and m, 
i.e. we report n/n 0 and m/n 0 . On the other hand, for 
the input parameters we set /3 as the main scale, i.e. we 
report the physical results as functions of 

Pn, /3/i, and A = y//3 g , (3) 

where the contact C is made dimensionless by Co, the 
unpolarized result at /3/r = 0. Note that, since we use 
an imaginary chemical potential difference in our Monte 
Carlo studies, the “bare” outcome of these simulations 
will be given as a function of /?/z> /3hi, and A. 

From the density and the magnetization one may ob¬ 
tain the isothermal compressibility and magnetic suscep¬ 
tibility simply by taking derivatives with respect to /3/r 
and /3/i, respectively. Mixed response functions (of n with 
respect to (3h, or m with respect to /3/i) may be obtained 
in the same fashion. On the other hand, numerical inte¬ 
gration of n with respect to /3/i provides the pressure for 
each value of fih. 
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The computational method utilized in this work is very 
similar to the one of Refs. |2T1 - [23| , but reduced to one spa¬ 
tial dimension and generalized to complex chemical po¬ 
tentials as explained above (see also Ref. [Ml)- Because 
one-dimensional problems are computationally inexpen¬ 
sive, it is possible to calculate in very large lattices, e.g. 
N x ~ 0(1O 2 ). For such sizes, the continuum limit is eas¬ 
ily achieved by lowering the density, while still remaining 
in the many-particle (i.e. thermodynamic) regime. For 
the proof-of-principle calculations presented here, we fix 
A = 1.0. This was chosen as being in the intermediate- 
to-strong-coupling regime, which is typically outside the 
range of validity of perturbative approaches. For such 
a coupling strength, a lattice size of N x = 61, which we 
fixed throughout this study, is sufficient to provide a good 
quantitative understanding of the continuum limit (see 
Ref. [24]). The physical extent of the system is L = N x i. 
where l = 1 fixes the spatial lattice units. The extent of 
the temporal lattice is given by /3 = tN t , where we take 
r = 0.05/£ 2 . 
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Figure 1. (Color online) Left: Density as a function of the 
imaginary chemical potential difference phi at various values 
of Pfi for a dimensionless coupling of A = 1.0. Right: Analytic 
continuation of the density as a function of Ph at various 
values of Pfi. In both plots the density is an even function 
about the origin and is plotted in units of the density of the 
noninteracting, unpolarized system. 


B. Computing the density and magnetization at 
imaginary asymmetry 


At imaginary chemical potential asymmetry, the par¬ 
tition function of Eq. ([2]) can be written in terms of a 
Hubbard-Stratonovich auxiliary held er as 

Z = y T>a\ det(l + zU[cr])\ 2 , (4) 


where z = = exp(/3 fiP) = exp(/3//T) = zf, and U[a\ is 

a matrix that encodes the dynamics of the system (see 
e.g. Ref. [TT]). We thus identify 

P[a\ = | det(l + zU[a})\ 2 (5) 


as the non-negative probability measure for our Monte 
Carlo calculations. The total (average) particle density 
n is then obtained in those calculations using 


1 5 In Z 

Ld(M 


— / VaP[a]Re 


Tr 


zU[o 


1 + zU[a\ J 


(6) 

and the (average) spin magnetization m can similarly be 
calculated using m = {\/L)d\w Z/d(Ph). To circumvent 
the sign problem in our Monte Carlo simulations, how¬ 
ever, we rather compute 


_ 1 din Z 
“ Ld(phi) 


2 i 


VaP[a]lm 


Tr 


( zU \°\ v 

\l + zU[a})_ ' 


(7) 

Since we assume that m is analytic as a function of gen¬ 
eral complex-valued h 7 at least in a finite domain about 
h = 0, we shall use the same label as for the physical spin 
magnetization. 

To determine the contact, we use the same approach 
as in Ref. [M]- The definition in ID is 


c _ 2 d<j30) 

PX T d(a 0 /X T ) 


h,t 


( 8 ) 


where D = — X l n Z is the grand thermodynamic poten¬ 
tial. Using the Feynman-Hellman theorem, it can be 
shown that 


C=-g(V), (9) 

where (V) is the thermal expectation value of the inter¬ 
action operator. The latter can be computed in Monte 
Carlo calculations using derivatives of ln Z with respect 
to the bare lattice coupling g or the lattice spacing r. 


III. RESULTS 

A. Monte Carlo results for imaginary asymmetry 

Throughout this work, we present Monte Carlo cal¬ 
culations at p = 8.0 and a lattice size of N x = 61. 
For the purposes of demonstrating the method, we fix 
the dimensionless coupling to A = 1.0, although a vari¬ 
ety of coupling strengths may also be explored using the 
same technique. The imaginary asymmetry parameter 
phi was varied over a full period [—7r, 7r], and the chem¬ 
ical potential parameter pfi was varied in the interval 
[—4.0,4.0], covering the semiclassical regime (where the 
virial expansion is valid) to the fully quantum mechanical 
regime. For each point in the plots below, we have taken 
1000 de-correlated Monte Carlo samples, thus ensuring 
that the statistical uncertainty is below 10%. 

In Fig. [l] we show the density as a function of phj and 
ph , respectively, for representative values of Pfi. Simi¬ 
larly, Fig. [2] displays the magnetization and Fig. [3] shows 
Tan’s contact. The statistical error of the Monte Carlo 
calculations on the imaginary side is estimated to be on 
the order of the symbol size in all three figures. In all 
cases we show the Monte Carlo results at Phi on the 
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Figure 2. (Color online) Left: Magnetization as a function 
of the imaginary chemical potential difference j3h\ at various 
values of /3/u for a dimensionless coupling of A = 1.0. Right: 
Analytic continuation of the magnetization as a function of 
j3h at various values of /3/x. In both plots the magnetization 
is an odd function about the origin and is plotted in units of 
the density of the noninteracting, unpolarized system. 
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Figure 3. (Color online) Left: Tan’s contact as a function 
of the imaginary chemical potential difference j3h\ at various 
values of (5fi for a dimensionless coupling of A = 1.0. C 0 is 
the contact at /3h = /3fi = 0. Right: Analytic continuation 
of the contact as a function of /3h at various values of /3/z. In 
both plots the contact is an even function about the origin. 
The curves on the right are color-wise paired by their value of 
P/j,, which coincides with the value on the left (by color code, 
or from top to bottom). Dashed-dotted lines give the results 
from a fit of the data to the polynomial-type ansatz ( |16| ), 
whereas solid lines result from fits of the data to our Pade- 
type ansatz (141 for the fit functions. 


left panel, and the corresponding analytic continuation 
(described in detail below) on the right. Although the 
imaginary side of the problem is not physically meaning¬ 
ful, the results are non-perturbative, and it is reassuring 
that the data falls on smooth curves that respect the even 
or odd symmetry around /3h\ = 0. We therefore display 
only the positive interval f5h\ £ [0 , tt] . 


B. Analytic continuation to the real axis 


In order to obtain the results of physical interest, we 
need to analytically continue the data from our Monte 
Carlo study to real-valued chemical potential differences. 
To this end, we fit our data to a specific ansatz for, e.g., 
the spin magnetization. Clearly, this is a critical step as 
the functional form of the ansatz is a priori unknown. 
However, as already discussed above, we know that, de¬ 
pending on the observable, the ansatz must be either an 
odd or even function in /3hi and that the partition func¬ 
tion is periodic in /3hj. Moreover, the virial expansion of 
the partition function Z suggests that Z can be written 
as a (asymptotic) series in powers of cos(/3/ii), see also 
our discussion in Sect. |IV| From such an analysis of the 
virial expansion, it also follows that higher-order terms 
in cos(/3/ji) become particularly important for fin 1. 
In order to take such higher order terms effectively into 
account, an ansatz of the type 


1 


1 + cos(/3hi) 


( 10 ) 


may be considered appropriate as it can be rewritten as 
an asymptotic series in powers of cos(/3/ii), which still 
obeys the periodicity in /3h\. An ansatz of this type may 
be viewed as a Pade approximant of cosines. 

To be specific, we choose to fit our Monte Carlo data 
for the magnetization at finite /3h\ with the function 


f(Ph i) = -i 


Msin(/3/ii) 


1 + B cos(/?/ii) + CF{f3h\) 


( 11 ) 


where A, B 1 and C are free real-valued fit parameters and 
the function F is given by 


F(x) = cos (x 3 /7r 2 ) . 


( 12 ) 


Note that this last function has been found empirically. 
As written, it is analytic on the whole complex plane, 
in particular on a disk of radius n centered at the ori¬ 
gin. However, it is not periodic in j3h\, as we required 
above. Here, we have indeed given up this constraint as 
we found that a large number of fit parameters is required 
to meet this criterion, rendering the fitting algorithm po¬ 
tentially unstable. For example, an ansatz of the form 
~ X)fc=i sm(fe/3/ii)/(l + J2k=i cos(k/3hi)) is com¬ 
patible with the above-mentioned constraints. However, 
the fit is of much lower quality than those obtained using 
Eq. Interestingly, the parameter C associated with 

the function F is in most cases found to be small when 


we fit our Monte Carlo data with the ansatz (|11|), see 
Tab. HD 


The analytic continuation f{/3h) of Eq. (11) is obtained 


simply by setting hi = — ih , and is therefore given by 


m = - 


Msinh(^Zi) 


1 + B cosh (/3/i) + CF{j3h) 


(13) 


where F is the analytic continuation of F. 
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In the case of the density and the contact, on the other 
hand, we expect even functions of fihi, and therefore we 
have chosen to fit the function 


g(Phi) = 7 


1 + Ar)(fih\) 


(14) 


where r/(x) = 1 — cos(a;). Here, 7 is a parameter that is 
fixed by the exactly known value at fih = fih\ = 0 , while 
A, B and C are free real-valued fit parameters. Once the 
parameters are obtained, the analytic continuation g(fih) 
of Eq. (141 is given by 


g(fih) = 7 


1 + Arj(fih) 


1 + Bij(fih) + Ci)((fih) 3 /ir 1 2 )_ 


(15) 


where ff{x ) = 1 — cosh(cc). 

The fit parameters for the density, magnetization, and 
Tan’s contact are provided in Tables |TJ|TTJ and 1 11 1[ respec¬ 
tively. Since the fits to the Monte Carlo data are sensitive 
to initial parameter values, the fitting algorithm performs 
several such fits with random initial parameter values in 
the interval [— 1 , 1 ], and the best fit with the minimum 
mean residuals is chosen for the analytic continuation. 
Given the functional form of Eqs. (Ill and (141, one 
should consider that poles may appear in either the fit 
or analytic continuation for a given set of fit parameters. 
Since we expect these quantities to be analytic for one di¬ 
mensional Fermi gasesjjthe fitting algorithm eliminates 
any fits that demonstrate such behavior. Of course other 
functional forms may be considered for the analytic con¬ 
tinuation with consideration to the constraints discussed. 

In Fig. [I] we also show the analytic continuation of the 
density, in Fig. [2] of the magnetization, and in Fig. [3] 
of Tan’s contact to real-valued chemical potential differ¬ 
ences. Although these functions are not periodic in fih, 
they remain valid only in the original restricted domain 
of [— 7 r, 7 r]. A few representative values of fin are shown 
in Figs. 00 and[3j however such analytic continuations 
may be performed for many values of fin on an unre¬ 
stricted domain, and the equations of state for various 
imaginary asymmetries may be constructed. Such plots 
for representative values of fih at a dimensionless cou¬ 
pling strength of A = 1.0 for the density, magnetization, 
and Tan’s contact are shown in Section HVBI 

One of the most interesting features we observe in our 
results is the behavior of the magnetization m/n 0 as a 
function of fih and fin- On the imaginary side (at least 
in the region studied), this quantity is non-monotonic 
in both of those variables. In particular, we note that 
the ordering of the curves, for different values of fin , is 


1 In higher dimensions, phase transitions may occur, potentially 

rendering physical quantities non-analytic at the transition point. 

Depending on the observable and the employed ansatz, such 
poles may therefore not just be an artifact resulting from the 
details of the fitting procedure but may be a hint to an underly¬ 
ing physical effect. 


Table I. Fit parameters for the density as they appear in 
Eq. (141 at a constant dimensionless coupling of A = 1.0 for 
various values of fin-: as well as the \ 2 per degree of freedom 
for each fit. Note that 7 is not a fit parameter (see main text). 
The value in parentheses indicates the calculated uncertainty 
of the least significant digit for each fit parameter. 


fin 

7 

A 

B 

C 

x 2 

-3.6 

1.0080 

-0.9670(4) 

-0.041(1) 

-0.001(1) 

0.71 

-3.2 

1.0172 

-0.9511(4) 

-0.0565(9) 

-0.002(1) 

0.58 

-2.8 

1.0294 

-0.9281(4) 

-0.080(1) 

-0.006(1) 

0.54 

-2.4 

1.0461 

-0.8942(9) 

-0.118(3) 

0.005(3) 

1.29 

-2.0 

1.0698 

-0.848(1) 

-0.156(3) 

-0.000(4) 

1.66 

-1.6 

1.1026 

-0.7904(6) 

-0.198(2) 

-0.016(2) 

0.35 

-1.2 

1.1245 

-0.716(1) 

-0.276(3) 

-0.008(4) 

1.25 

-0.8 

1.1731 

-0.6334(6) 

-0.323(2) 

-0.019(3) 

0.63 

-0.4 

1.1921 

-0.5433(9) 

-0.364(4) 

-0.034(7) 

1.30 

0.0 

1.2240 

-0.492(2) 

-0.415(2) 

-0.056(5) 

2.74 

0.4 

1.2359 

-0.503(3) 

-0.493(3) 

-0.009(1) 

5.21 

0.8 

1.2361 

-0.4(1) 

-0.4(1) 

0.00(1) 

0.54 

1.2 

1.2214 

-0.18(5) 

-0.22(5) 

-0.013(6) 

0.97 

1.6 

1.2109 

-0.32(2) 

-0.34(2) 

0.004(1) 

0.71 

2.0 

1.1948 

-0.33(3) 

-0.35(3) 

0.005(2) 

1.81 

2.4 

1.1763 

-0.22(2) 

-0.24(2) 

-0.002(1) 

0.41 

2.8 

1.1625 

-0.11(3) 

-0.13(3) 

-0.007(1) 

0.43 

3.2 

1.1507 

-0.20(3) 

-0.22(3) 

-0.001(1) 

0.45 

3.6 

1.1397 

-0.07(6) 

-0.08(6) 

-0.009(3) 

0.60 


partially inverted at large enough fih. This behavior, 
however, results in a perfectly ordered set of curves on 
the real-valued (fih) side, in a way that respects both 
thermodynamic stability and physical intuition. 

In Fig. [3] we show two possible fits and their corre¬ 
sponding analytic continuations for the contact, namely 
Eqs. ( |T4| ), ( fl5| ), and an alternative function 

q(fihi) = 7 [l + An(fihi) + Br]((fihi) 3 /Tr 2 )\ , (16) 


where again A and B are free real-valued fit parameters 
and 7 is a fixed value, as discussed previously. The ana¬ 
lytic continuation q(fih) is given in terms of fj(fih). While 
the fits on the imaginary side appear to be of compara¬ 
ble quality, they differ enough in the details that their 
analytic continuation to the real side displays noticeable 
discrepancies. This is particularly evident for large fin- 
We take this to be indicative of the limitations of our 
approach and it should be viewed as a warning with re¬ 
spect to the choice of the fit function: The Pade form 
given in Eq. (14) effectively takes into account arbitrarily 
high powers of cos(fih \), whereas the ansatz (16) may be 
viewed as a low-order approximation of the ansatz (14) 
in powers of cos(fih\) which is expected to be valid only 
in the vicinity of fih\ = 0. A priori , it is difficult to 
judge under which conditions a low-order approximation 
is justified at all. In the present case, for example, the 
value of the coupling does not provide a direct criterion. 
In fact, already the free Fermi gas in one dimension for 
fin > 0 is described by an asymptotic series in powers 
of cos(fihi). For fin <C —1, on the other hand, it can 
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Table II. Fit parameters for the magnetization as they appear Table III. Fit parameters for the contact as they appear in 


in Eq. {Ill at a constant dimensionless coupling of A = 1.0 for Eq. (141 at a constant dimensionless coupling of A = 1.0 for 


various values of /3/x, as well as the y per degree of freedom 
for each fit. The value in parentheses indicates the calculated 
uncertainty of the least significant digit for each fit parameter. 


various values of /3/r, as well as the y per degree of freedom 
for each fit. Note that 7 is not a fit parameter (see main text). 
The value in parentheses indicates the calculated uncertainty 



A 

B 

C 

x J 

of the least significant digit for each fit parameter. 


-3.6 

1.022(4) 

0.044(3) 

0.011(5) 

0.48 


7 

A 

B 

C 

x 2 

-3.2 

1.031(4) 

0.062(3) 

0.011(5) 

0.50 

-3.6 

0.0003 

-0.4(2) 

-0.60(4) 

0.15(4) 

0.47 

-2.8 

1.048(3) 

0.094(3) 

0.012(4) 

0.41 

-3.2 

0.0019 

-0.50(5) 

-0.50(3) 

-0.00(5) 

0.52 

-2.4 

1.047(6) 

0.146(5) 

-0.015(6) 

1.01 

-2.8 

0.0062 

-0.503(2) 

-0.10(5) 

-0.40(8) 

0.12 

-2.0 

1.062(5) 

0.202(4) 

-0.023(6) 

1.35 

-2.4 

0.0159 

-0.5(1) 

-0.40(9) 

0.4(6) 

1.01 

-1.6 

1.098(4) 

0.269(3) 

0.000(4) 

0.48 

-2.0 

0.0390 

-0.45(8) 

-0.39(6) 

0.3(4) 

1.30 

-1.2 

1.134(4) 

0.387(3) 

0.017(4) 

1.80 

-1.6 

0.0971 

-0.5025(5) 

-0.22(2) 

-0.27(3) 

0.42 

-0.8 

1.142(3) 

0.521(3) 

0.037(3) 

5.81 

-1.2 

0.1712 

-0.49(1) 

-0.42(1) 

-0.00(4) 

0.98 

-0.4 

1.089(4) 

0.632(3) 

0.052(3) 

8.90 

-0.8 

0.3542 

-0.496(2) 

-0.418(4) 

-0.05(1) 

0.32 

0.0 

0.978(3) 

0.702(2) 

0.067(2) 

23.32 

-0.4 

0.5722 

-0.493(2) 

-0.473(3) 

0.006(8) 

0.64 

0.4 

0.833(3) 

0.711(3) 

0.088(3) 

36.24 

0.0 

1.0000 

-0.496(1) 

-0.464(3) 

-0.016(6) 

1.02 

0.8 

0.695(3) 

0.684(3) 

0.123(3) 

7.36 

0.4 

1.4929 

-0.4974(7) 

-0.450(2) 

-0.029(4) 

0.70 

1.2 

0.574(2) 

0.660(3) 

0.156(3) 

4.95 

0.8 

2.1795 

-0.4978(7) 

-0.402(2) 

-0.075(4) 

0.60 

1.6 

0.473(2) 

0.651(3) 

0.175(3) 

7.57 

1.2 

2.7994 

-0.494(1) 

-0.386(2) 

-0.061(6) 

0.49 

2.0 

0.397(2) 

0.639(4) 

0.193(4) 

8.02 

1.6 

3.6330 

-0.4986(7) 

-0.331(3) 

-0.138(7) 

1.09 

2.4 

0.340(2) 

0.630(5) 

0.214(5) 

11.10 

2.0 

4.3998 

-0.498(1) 

-0.304(4) 

-0.16(1) 

2.21 

2.8 

0.292(2) 

0.634(6) 

0.215(6) 

14.55 

2.4 

5.0691 

-0.4980(5) 

-0.298(2) 

-0.157(5) 

0.61 

3.2 

0.256(1) 

0.638(5) 

0.218(4) 

14.73 

2.8 

5.8386 

-0.4975(5) 

-0.278(2) 

-0.166(5) 

0.53 

3.6 

0.227(2) 

0.636(7) 

0.218(7) 

24.89 

3.2 

6.6344 

-0.4967(7) 

-0.256(2) 

-0.180(7) 

0.69 






3.6 

7.3683 

-0.4982(6) 

-0.240(2) 

-0.211(6) 

0.85 


be shown that already a low-order approximation yields 
reliable results for the free Fermi gas, see also our dis¬ 
cussion of the virial expansion in Sect. IV B| Appar¬ 
ently, the strength of the coupling does not enter these 
arguments. A variation of the strength of the coupling 
is only expected to change the numerical values of the 
series coefficients associated with such an expansion in 
powers of cos(/3hi) and may therefore only effectively im¬ 
prove or worsen the convergence properties of this series. 
Note that both limits /?/i 1 and /3fi <C —1 correspond 

to weak-coupling limits in our Monte Carlo study with 
fixed dimensionless coupling A and fixed inverse temper¬ 
ature /3. In the regime |/3^t| < 1, on the other hand, the 
theory is effectively in the strongly coupled regime, see 
also our discussion Sect. IIV A1 


C. Magnetization-to-density ratio and magnetic 
susceptibility 

It is important to understand whether the system we 
are studying is appreciably magnetized in the region of 
parameter space that we explore here. To clarify this 
point, in Fig. [4] we show the ratio of the magnetization 
m to the density n. In absolute value, this ratio can 
only vary between 0 (unpolarized) and 1 (fully polar¬ 
ized). Furthermore, it is reassuring that m/n lies within 
the interval [0,1] after analytic continuation, and is a 
monotonically increasing function with f3h. 

Our results for the magnetization allow us to compute 
the magnetic susceptibility y of a polarized Fermi gas 



Figure 4. (Color online) Ratio of the magnetization m to 
the density n as a function of real-valued /3h and /3fi = - 
2.0, -1.0, 0, 1.0, 2.0. The solid lines show the second-order 
virial expansion at pfx = —2 (top) and -1 (bottom). Note 
that the virial expansion works well for /3/x = — 2, but fails 
dramatically for pfj, = — 1 and above; see Section |IV B| for 
details. 


simply by taking a derivative: 


1 dm 
n 0 d(/3h) ’ 


(17) 


where in practice we simply take an analytic derivative 
of Eq. |I3| ) for each discrete value of /3/r. The magnetic 
susceptibility as a function of /3h for representative val¬ 
ues of /3fi at a fixed dimensionless coupling of A = 1.0, 
as well as the noninteracting case is shown in Fig. [5] 
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Figure 5. (Color online) The magnetic susceptibility \ as 
a function of /3h at representative values of /3p, obtained by 
taking an analytic derivative of the polarization with respect 
to j3h. Left panel: interacting case at A = 1. Right panel: 
noninteracting case in the continuum. 


Note that the deterioration in the accuracy of the ana¬ 
lytic continuation at large (3h (as discussed in Sec. Ill B) 
becomes more apparent as we take derivatives of physical 
observables. 



ph i j8/i 


IV. COMPARISON WITH OTHER 
APPROACHES 

In this section, we compare the results from our Monte 
Carlo simulations with those from other approaches 
which also helped us to guide the analytic continuation 
of our data from imaginary to real-valued chemical po¬ 
tential differences. Moreover, these comparisons allow 
us to gain at least some insight into the phenomenology 
underlying one-dimensional Fermi gases. 


A. Pairing effects 

The partition function of the noninteracting gas (A = 
0 ) may be computed analytically in the grand-canonical 
ensemble using 


Figure 6. (Color online) Top: Density n of the noninteracting 
Fermi gas for several values of /3p as a function of imaginary 
(left panel) and real (right panel) values of /3h. Bottom: Same 
as top panel, but showing the magnetization m. Both n and m 
are displayed in units of the density n 0 of the noninteracting, 
unpolarized system. 


which, in retrospect, motivates our general forms for the 
ansatze for the fit functions used to analytically con¬ 
tinue the Monte Carlo data from imaginary to real-valued 
chemical potential differences. The coefficients bk can be 
related to the one-, two-, three-, ..., IV-body problem, 
see also our discussion below. Note that this series con¬ 
verges particularly well for /3/z <C —1. 

The density and spin magnetization follows immedi¬ 
ately from Eq. (18). For the density, we obtain 


lnZ(/3 n,fih) = —t— [40 t ) + 4( 2 J.)] . (18) 

Y 7T Aji 

where 

/ OO . . 

dxhi M + z s e~ x2 J , (19) 

-OO ' ' 

and and z ^ = e^e - ^. From this, it fol¬ 

lows immediately that In Z can be written in terms of an 
asymptotic series of the form 

OO 

In Z(pn, 0h) = cosh (kph ), (20) 

k=0 


x A t <91nZ(/3/i, fih) 

= T d(M = 

where I^Zg) = z s dl 0 (z s )/dz s , 
tion m, we find 


A= i/i(*t)+ j i (*t)]. 
v 71 " 

( 21 ) 

and, for the magnetiza- 


. A t d\nZ(/3fi, fih) 

mXT = T am 



( 22 ) 


In Fig. [6j we show our results for the density and the mag¬ 
netization of the free Fermi gas as a function of $h\ for 
various values of fifi. Comparing these results for the free 
Fermi gas with those from our Monte Carlo study (see 
Figs, [l] and [2]), we observe that both agree qualitatively, 




























at least for finite (3p. For (3p = 0, on the other hand, 
we observe that the results from the free gas diverge for 
\/3hi\ —> 7 r. This can be understood from Eq. ( fl9| ): Set¬ 
ting f3p = 0, we observe that the function I 0 diverges at 
least logarithmically in the limit |/3/q| —> tt. We would 
like to add that the reason for the divergence appear¬ 
ing in the results in this limit can also be understood 
from a study of the free propagator which becomes sin¬ 
gular for (3p = 0 and p 2 = 0 in the limit \f3hi\ —> n. 
Formulated in the language of thermal field theory, the 
fermionic Matsubara modes /3w n = (2n+l)7r entering the 
free propagator effectively assume the form (3ui n = 2m: 
associated with bosonic degrees of freedom in the limit 
|/3 hi | —» 7r, see also Ref. [13]. In other words, in this limit 
the fermions acquire a (thermal) zero mode. However, 
we emphasize that the partition function is analytic for 
any finite value of f3pr\ Therefore, no divergences oc¬ 
cur in the limit \/3hi\ —> 7r for f3p ^ 0. The situation is 
substantially different in the case of an interacting Fermi 
gas as studied with our Monte Carlo approach. Here, (3p 
is only a parameter with limited physical meaning. In 
fact, whereas /i determines the energy of the free Fermi 
gas, the parameter p in the Monte Carlo study only sets 
the scale for the energy of the corresponding interacting 
Fermi gas. Loosely speaking, an effective chemical po¬ 
tential p e g may also be assigned to the interacting Fermi 
gas. In general, its value would then be different from 
the value of the parameter p. Indeed, even for f3p = 0, 
the results for physical observables from our Monte Carlo 
study appear to be analytic as a function of (3h\. For ex¬ 
ample, the spin magnetization m diverges for f3p = 0 
and \(3hi\ — > tt in the case of the free Fermi gas, see 
Fig. [6] For the interacting Fermi gas, on the other hand, 
to appears to be analytic for f3p = 0 and only exhibits 
a rapid decrease in the limit \f3h\\ -4 7r, suggesting the 
effective chemical potential p e g associated with the in¬ 
teracting theory is small but finite, see Fig. [5] 

Although there is no spontaneous symmetry break¬ 
ing in one dimension in the long-range limit, pairing 
of fermions is a priori still possible for any finite cou¬ 
pling strength and expected to impact the ground state. 
It is therefore worthwhile to study pairing in the one¬ 
dimensional Fermi gas at zero temperature. In this case, 
the coupling is conveniently rendered dimensionless with 
the aid of the chemical potential p which sets the scale: 

9 = = -4=f ■ (23) 

vWl V\iA 

Thus, the scale /3 drops out as it should be and the di¬ 
mensionless coupling g can now be viewed as a measure of 
the potential energy (measured in terms of g) relative to 
the kinetic energy (measured in terms of p). We observe 
that, for fixed A and /3 as in our case, we approach the 


2 Note that In Z is a sum of /ofe^e^) and /ofe^e & h ) up to 
numerical factors. 



Figure 7. (Color online) Dimensionless binding energy Eb/ p 
of the two-body bound state in the presence of (inert) Fermi 
surfaces as a function of h/p for g = tt. The gray-shaded area 
depicts the regime in which the formation of a bound state 
with finite center-of-mass momentum is favored. 


weak-coupling limit for, e.g., f3p 1. On the other hand, 
the theory becomes strongly coupled for fixed A ~ 0(1) 

if |/3/t| < 1. 

Keeping this in mind, let us now analyze the role of 
pairing effects in our Monte Carlo study by simply con¬ 
sidering the two-body problem in the presence of two 
Fermi surfaces, in close analogy to standard BCS the¬ 
ory [25] . The underlying Schrodinger equation, which 
has proven very useful to understand the general phase 
structure of imbalanced Fermi gases nnna, is given by 


[ e s( d x s )-gS(x t -x l ) + E b T(aJ t ,X 4 .) = 0. (24) 


s=t4 


Here, T is the wave-function of the bound state. The 
operator e s is defined as e s (d x ) = \ — (2m)~ 1 d 2 — eF, s | 
with ep jS (s =t,4.) being the Fermi energy of the up- and 
down-fermions, respectively. Interestingly, the solution 
of this one-dimensional two-body problem can in princi¬ 
ple also be given in closed form m For our purposes, 
however, only the (binding) energy of the lowest-lying 
bound-state, which is obtained from a minimization of 
the energy Eb with respect to the total momentum P, is 
of particular interest. For illustration purposes, we show 
the energy of this state as a function of h/p in Fig. [7] 
for g = 7 r (in the strong-coupling regime), correspond¬ 
ing to [3p = 1 / 7 T 2 « 0.1 in our Monte-Carlo study with 
fixed A = 1.0. The gray-shaded area in Fig. [7] depicts the 
regime in which it is energetically most favorable to form 
a bound state with finite center-of-mass momentum. In 
a full many-body treatment, the true ground state can 
potentially be inhomogeneous in this regime [251 EH- 
We observe that the dimensionless binding energy 
Eb/ p becomes smaller for increasing spin-imbalance h/p. 
Moreover, we find that the formation of a two-body 
bound state is no longer energetically favored for g < 0.96 
which corresponds to (3p > 1.08 in our Monte Carlo study 
with fixed A = 1.0. Loosely speaking, this suggests that, 
for fixed A, the Fermi gas undergoes a crossover from a 
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strongly correlated to a weakly correlated system in the 
limit fig 1. 

A word of caution needs to be added here: Our study 
of bound-state formation in the presence of (inert) Fermi 
surfaces is clearly only an approximation as the Fermi 
surfaces are smeared out at finite temperature and cou¬ 
pling strength. Moreover, it has been found that the spin- 
balanced 1V| + A^-problem “dimerizes”, i.e. the ground 
state energy of this system in the zero-temperature limit 
is given by the (N+ + Nfi)/2 times the binding en¬ 
ergy of the associated two-body bound state (see, e.g., 
Refs. [3 I2HH33)- Thus, even in the limit of small di¬ 
mensionless coupling g , bound states are formed. Never¬ 
theless, as in our study of bound-state formation in the 
presence of Fermi surfaces, the dimensionless energy of 
the system (i.e. energy measured in units of g) decreases 
for fixed coupling g and increasing g. However, the crit¬ 
ical coupling for bound-state formation turns out to be 
zero in the exact solution, independent of the degree of 
spin imbalance of the system. 

With respect to our Monte Carlo simulations, these 
considerations imply that pairing effects are present for 
all values of g and h and are at least partly responsible for 
the difference between the results for the free Fermi gas 
and our Monte Carlo results. It should also be added that 
our Monte Carlo study is bound to finite temperature 
\fih\ < 7T which hinders a direct comparison between our 
Monte Carlo results and the exact solutions [3 HE] only 
available for the zero-temperature limit. At finite tem¬ 
perature, thermal energy is “pumped” into the system, 
effectively resulting in dissociation of the bound states^] 

We close by adding a comment on local ordering in 
one-dimensional Fermi gases. The formation of bound 
states can be considered as a necessary condition for 
the formation of a superfluid condensate. Of course, as 
stated above, there is no spontaneous symmetry break¬ 
ing in these one-dimensional systems in the long-range 
limit. Nevertheless, the emergence of local ordering, i.e. 
the emergence of a condensate in the presence of a suf¬ 
ficiently large infrared cutoff, may be possible. In three- 
dimensional systems, these types of phases are associated 
with precondensation EH ED- In experiments, such an 
infrared cutoff scale is effectively set by the inverse of the 
length scale associated with the confining geometry, e.g. 
a harmonic trap potential. Since the extent of inhomoge¬ 
neous phases in the space spanned by the experimental 
parameters is expected to be large [251E2, it may indeed 
be worthwhile to further study the fate of these phases 
at finite temperature with the aid of our present Monte 
Carlo setup. 


3 Note that our Monte Carlo results are always given as a function 
of 0h. Thus, large values of 0h correspond to low temperatures 
for fixed chemical potential difference h, whereas small 0h is 
associated with high temperatures in this case. 


B. Virial expansion 

Let us now consider the virial expansion of the parti¬ 
tion function, i.e. an expansion in powers of z = e^. In 
the fig —> —oo limit, where the virial expansion is valid, 
we can evaluate the density and the magnetization order 
by order. Indeed, at leading order in z, n^A T = z^ 
and therefore 

n\ T = (n^ + n^)\ T = 2e^ M cosh(/3/i) (25) 
which leads to 

77 

— = cosh(/3/i) (26) 

n 0 

where n 0 is the density for the unpolarized system. This 
result is the leading order in the virial expansion and does 
not depend on the interaction. Similarly, we find for the 
magnetization that 

777, 77- a . — 71 1 

— = -1 -± = sinh(fih) , (27) 

n o n 0 

at leading order in z which then yields 

777 

— = tanh (fih) , (28) 

n 

which is valid at the same (leading) order in z. 

In general, accessing higher orders in the virial expan¬ 
sion requires solving the two-, three-, ..., ./V-body prob¬ 
lems (see, e.g., Ref. 33] ). The grand-canonical partition 
function for systems with chemical potential asymmetry 
may be written as 


Z=J2Y1 zNwM QN t ,N j. , 


(29) 


N =0 M 


where N = is the total particle number and 

M = measures the spin polarization. Moreover, 

we have introduced the quatitites w = e^ h and Qn t ,n 1 
being the canonical partition function of a system with 
Ah spin-up fermions and N, spin-down fermions. Ex¬ 


panding Eq. (29) to second order in z yields 


Z = 1 + 2Q X 0 z cosh(fih) 

+ Q ia z 2 + 2Q 2 ,o~ 2 cosh(2 fih) + 0(z 3 ), (30) 

where Qyo = L/Xt and Q±g and Q 2,0 may be deter¬ 
mined by direct calculation [22] . Note that the expansion 
coefficients Q /v t .vy depend implicitly on the coupling A. 
The density n and magnetization m in the second-order 
virial expansion may then be determined from Z through 
Eqs. (211 and (22) and compared with Monte Carlo re¬ 
sults in the z —>■ 0 limit. We do this explicitly in Fig. 
and find excellent agreement for small z and fih; however, 
the quality of the agreement deteriorates quickly as fih 
is increased unless z ~ 0, as expected. For completeness, 
Fig.[9]shows Tan’s contact as a function of fig, for several 
values of fih , and at A = 1. 
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Figure 9. (Color online) Analytic continuation of Tan’s 
contact as a function of /3/j, at a constant coupling A = 1.0 
and /3h = 0, 0.5, 1.0, 1.5, 2.0. The statistical error is on the 
order of the size of the plotted symbol. C n is the contact at 
/3n = /3h = 0. 

of the free gas in Fig. [lO] In all cases the results ob¬ 
tained with our proposed analytic-continuation approach 
lie generally between the leading (noninteracting) calcu¬ 
lation and the next-to-leading-order result. 


Figure 8. (Color online) Top: Analytic continuation of the 
density n (in units of its unpolarized, noninteracting counter¬ 
part n 0 ) as a function of /?/* at a constant coupling A = 1.0 
and /3h = 0, 0.5, 1.0, 1.5, 2.0. The solid black line shows the 
second-order virial expansion for each value of /3h. Bottom: 
Same as top, but showing the magnetization m. Error bars 
were estimated by varying the fit parameters by an amount 
given by the uncertainty in the calculated fits. 


C. Lattice perturbation theory 


As an additional verification of the equations of state 
obtained through an analytic continuation, we performed 
a next-to-leading order lattice perturbation theory cal¬ 
culation by expanding the lattice grand-canonical parti¬ 
tion function Z(/3n,/3h) in the dimensionless parameter 
A = 2{e T9 - 1) (which arises naturally in lattice cal¬ 

culations, see Ref. El) about the noninteracting limit 
(A = 0), up to the second non-vanishing term. Such an 
expansion on the lattice yields 


In Z = In Z( 



where Zo(f3fjb, (3h) is the partition function of the nonin¬ 
teracting gas, e k = k 2 /2m and the sum over k is over all 
possible lattice momenta. The density n/n o and magne¬ 
tization m/no hr terms of this perturbation theory there¬ 
fore follow from Eq. (31) using Eqs. (21) and ( [22] ). We 
display the results of this analytic calculation with the 
numerical Monte Carlo results and the equations of state 


V. SUMMARY AND CONCLUSIONS 

In summary, we have performed a non-perturbative 
characterization of the density n, magnetization m, mag¬ 
netic susceptibility y, and Tan’s contact C of a ID, at¬ 
tractively interacting Fermi gas. To this end, we imple¬ 
mented the conventional finite-temperature lattice Monte 
Carlo formalism, but generalized to include complex 
chemical potentials. When the chemical potential asym¬ 
metry h is purely imaginary, there is no sign problem 
and the Monte Carlo calculation can be carried out as 
usual. Our Monte Carlo results on the imaginary-/i side 
are therefore exact up to controlled (statistical and sys¬ 
tematic) uncertainties. 

To obtain results for real h, we performed fits to our 
numerical data and implemented an analytic continua¬ 
tion, i.e. we set h —¥ ih. In some regions of param¬ 
eter space, different functional forms for the fits may 
yield very different analytic continuations. However, very 
simple functional forms such as polynomials can be dis¬ 
carded as they are likely too simple to capture the essen¬ 
tial physics. Generally speaking, at low enough fih all 
fits lead to (approximately) the same analytic continua¬ 
tion. As we show in Fig. [4j however, even for low f3h, 
we achieve non-trivial magnetization ratios as large as 
m/n ~ 0.3 — 0.5. 

We have presented our results as a function of the di¬ 
mensionless parameters (3^ and /3/i, but focused on an 
intermediate- to strong-coupling regime A = \//3g ~ 0(1) 
as a non-trivial case of relevance for future studies. Our 
results for n and m agree qualitatively with the free 
Fermi gas where differences may be traced back to pair- 
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Figure 10. (Color online) Comparison of the Monte Carlo 
results after the analytic continuation, a next-to-leading order 
perturbation theory calculation (NLO PT), and the free gas 
for the density (top) and the magnetization (bottom). Results 
are displayed for /3h = 0.0,0.5,1.0,1.5 and 2.0 (bottom to top) 
in the strongly interacting region > —1. 
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Figure 11. (Color online) Analytic continutation of the 
density (top) and the magnetization (bottom) for the non¬ 
interacting system as a function of /3/r for various values of 
j3h. The exact solutions in the continuum limit are shown as 
solid black lines. 


ing effects. Moreover, our results are consistent with the 
second-order virial expansion, which is non-perturbative 
in the interaction, in the regime /3 n < 0, where that ex¬ 
pansion is valid. We also note, however, that the virial 
expansion deteriorates as flh is increased (at fixed z). 

The calculations carried out in this work correspond 
to fixed lattice volume L = N x £ = 61 and extent of the 
imaginary-time direction /3 = N t t = 8.0, where i = \ 
and r = 0.05. The associated systematic effects should, 
in principle, be further investigated, although the results 
of our previous work [23] indicate that those effects are 
below 10%. We consider our present work as a proof- 
of-principle study of our imaginary spin-imbalance ap¬ 
proach. However, results for, e.g., the thermal equations 
of state, could already be extracted from it and com¬ 
pared to present and future experiments IHj , if avail¬ 
able. In any case, our present study is mostly aimed 
at paving the way to more computationally demanding 
systems in two and three dimensions. For example, our 
approach allows one to map out to some extent the finite- 
temperature phase diagram of spin-imbalanced unitary 
Fermi gases in three dimensions, and therefore permits 
one to, at least, narrow down the regime in parameter 
space in which the critical point is located. In that re¬ 


gard, all of our present results indicate that calculations 
in higher dimensions should be feasible with the proposed 
method. 
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Table IV. Fit parameters for the density at alternative cou¬ 
pling strengths for various values of /3/r., as well as the X' 2 P er 
degree of freedom for each fit. Note that 7 is not a fit param¬ 
eter (see main text). The value in parentheses indicates the 
calculated uncertainty of the least significant digit for each fit 
parameter. 


A 

Pn 

7 

A 

B 

C 

x J 

0.5 

-2.0 

1.0272 

-0.8876(3) 

-0.1553(7) 

-0.0029(7) 

1.45 

0.5 

-1.0 

1.0624 

-0.7412(8) 

-0.299(2) 

-0.027(2) 

2.63 

0.5 

0.0 

1.0937 

-0.5216(4) 

-0.392(3) 

-0.064(4) 

2.21 

0.5 

1.0 

1.1018 

-0.2(3) 

-0.2(3) 

-0.05(6) 

0.73 

0.5 

2.0 

1.0908 

-0.20(2) 

-0.23(2) 

-0.014(2) 

0.64 

2.0 

-2.0 

1.2530 

-0.712(2) 

-0.163(7) 

-0.007(9) 

0.69 

2.0 

-1.0 

1.5051 

-0.483(6) 

-0.278(8) 

0.11(6) 

1.97 

2.0 

0.0 

1.6588 

-0.502(2) 

-0.491(2) 

-0.012(1) 

0.98 

2.0 

1.0 

1.5738 

-0.35(5) 

-0.37(5) 

0.006(1) 

1.78 

2.0 

2.0 

1.4263 

-0.17(5) 

-0.18(5) 

0.005(1) 

1.14 


Table V. Fit parameters for the magnetization at alternative 
coupling strengths for various values of p/j,, as well as the y 2 
per degree of freedom for each fit. The value in parentheses 
indicates the calculated uncertainty of the least significant 
digit for each fit parameter. 


A 

Pn 

A 

B 

C 

x 2 

0.5 

-2.0 

1.082(2) 

0.193(2) 

0.000(3) 

0.72 

0.5 

-1.0 

1.170(3) 

0.464(3) 

0.020(3) 

3.84 

0.5 

0.0 

1.109(3) 

0.778(2) 

0.094(2) 

26.24 

0.5 

1.0 

0.770(9) 

0.708(9) 

0.186(8) 

756.92 

0.5 

2.0 

0.487(7) 

0.64(1) 

0.25(1) 

458.75 

2.0 

-2.0 

1.069(5) 

0.179(4) 

0.010(6) 

0.51 

2.0 

-1.0 

1.023(8) 

0.392(7) 

0.027(8) 

1.57 

2.0 

0.0 

0.719(4) 

0.539(4) 

0.049(5) 

6.39 

2.0 

1.0 

0.422(2) 

0.562(3) 

0.093(4) 

9.16 

2.0 

2.0 

0.265(1) 

0.586(3) 

0.115(3) 

16.25 


eters A, B, and C are not necessarily smooth functions of 
phi, and noise is introduced into the result for m{Pn)/no. 
The smoothness of the interacting results shown in the 
main text for A > 0 significantly reduces such effects. 


Appendix B: Equations of state for alternative 
coupling strengths 


In the main text we have illustrated the method of 
calculating the equation of state of polarized, interacting 
fermions using complex chemical potentials at a constant 
coupling of A = 1.0 for clarity and as a proof of princi¬ 
ple. This same method can be applied to other values of 
the coupling as well, considering that larger values of A 
may require a larger number of Monte Carlo samples in 
order to mantain the statistical quality of the data. In 
Figs. [l2l and [l3| we show the density and magnetization 
as functions of both /3 h and ph\ for dimensionless cou¬ 
plings A = 0.5 and 2.0. The fit parameters for a chosen 
set of pfi at these couplings are provided in Tables IV 
andlVl 



0 0.5 I 1.5 
p/ h 


2 2.5 3 


0 0.5 1 


1.5 

jS/i 


2 2.5 3 


Appendix A: Analytic continuation of the 
non-interacting system 

In order to verify that the approach of analytically 
continuing the equations of state of the density n/n.Q 
and magnetization m/no is valid using the fit ansatze 
developed in the main text, we have performed the same 
prescription for the non-interacting polarized Fermi gas 
and compared with the exact solution for these quanti¬ 
ties in the continuum limit. The results comparing the 
Monte Carlo and analytic solutions are shown in Fig. 

There is excellent agreement between the two cal¬ 
culations within systematic and statistical error, which 
demonstrates a measure of validity for the analytic con- 
tinution in the interacting case. In the limit of A —> 0 on 
the imaginary side, m{ph\)/no converges to a sawtooth 
wave of periodicity 27r for large positive /3/x, whose be¬ 
havior in the vicinity of phi = ±7r is difficult to capture 
for the given ansatz. As such, the determined fit param- 



Ph p/i 


Figure 12. (Color online) Top: Density at a decreased cou¬ 
pling strength of A = 0.5 as a function of the imaginary chemi¬ 
cal potential difference Phj (left) and its analytic continuation 
to ph (right). Bottom: Magnetization at a decreased coupling 
strength of A = 0.5 as a function of the imaginary chemical 
potential difference Phi (left) and its analytic continuation to 
ph (right). 
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Vh pi> 


Figure 13. (Color online) Top: Density at an increased cou¬ 
pling strength of A = 2.0 as a function of the imaginary chemi¬ 
cal potential difference j3h\ (left) and its analytic continuation 
to /3h (right). Bottom: Magnetization at a increased coupling 
strength of A = 2.0 as a function of the imaginary chemical 
potential difference /3hi (left) and its analytic continuation to 
f)h (right). 
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